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SECTION  1 

TECHNICAL  D  SCUsSION 


While  running  the  S-CUBED  final  prediction  for  MISTY  ECHO,  it  was  noted 
that  the  on-axis  profiles  appeared  to  be  ringing  at  and  behind  the  front.  This 
manifested  itself  in  the  monitor  plots,  as  well,  v/ith  noticeable  ringing.  Following  the 
conclusion  of  the  prediction  calculation,  we  undertook  a  numerical  study  aimed  at 
improving  the  damping  in  STREAK.  This  report  documents  that  study. 

The  need  for  a  certain  amount  of  damping  in  finite  difference  solutions  of 
stress  wave  propagation  has  long  been  recognized.  It  is  also  generally  understood  that 
this  damping  may  be  explicit,  as  in  the  Richtmeyer-von  Neumann  artificial  viscosity,  or 
inherent  in  the  finite  difference  approximations  as  truncation  error,  as  is  the  case  in 
the  Eulerian  first  order  donor  scheme. 

The  need  for  damping  (code  glue,  to  some)  is  somewhat  disturbing,  but  can  be 
shown  to  have  both  physical  and  mathematical  justifications.  The  case  of  shock 
propagation  is  where  the  need  for  damping  is  most  clear.  Since  real  viscous  effects 
are  important  for  the  structure  of  physical  shocks,  it  should  not  be  surprising  that  a 
numerical  analogue  of  viscosity  is  required  for  numerical  simulation  of  a  shock.  There 
is,  however,  a  much  more  profound  explanation.  If  we  assume  a  steady  wave  with  a 
constant  profile,  it  can  be  shown  that  the  physical  loading  path  of  a  Lagrangian 
particle  is  along  a  Rayleigh  line.  (This  is  not  widely  appreciated,  but  it  is  true.)  The 
terminal  state  of  this  loading  is  the  Rankine-Hugoniot  state,  an  equilibrium  state,  but 
other  states  on  the  loading  path  are  non-equilibrium  states  which  are  not  compatible 
with  the  equation  of  state.  In  other  words,  the  material  states  required  by  the 
dynamics  of  shock  wave  propagation  are  not  consistent  with  the  equation  of  state.  In 
addition,  unless  there  are  non-equilibrium  contributions  to  the  pressure  which  can 
make  the  material  pressure  compatible  with  the  dynamics,  there  will  not  be  a 
mathematical  solution  to  the  steady  wave  problem.  This,  we  believe,  is  the  cause  of 
numerical  oscillation. 
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Thus  in  order  to  have  well-behaved  numerical  solutions,  it  is  essential  to  have 
viscous  or  other  non-equilibrium  stresses.  Low  order  schemes  produce  these 
contributions  by  accident.  High  order  schemes  require  an  explicit  treatment. 

In  the  early  1980’s  STREAK  was  upgraded  from  a  first  order  code  with  no 
explicit  damping  to  a  second  order  code  in  space  and  time  with  a  FRAM-like 
damping.^  This  scheme  looks  for  local  oscillations  and  introduces  damping  as  needed 
only  to  maintain  monotonicity.  The  code  was  extensively  tested  and  tie  FRAM 
method  was  found  to  perform  well  in  a  large  number  of  test  problems.  However, 
recent  production  calculations  with  STREAK  have  exhibited  undesirable  spikes  and 
ringing  which  have  been  traced  to  insufficient  damping  at  the  wavefront  when  the 
overall  computational  time  step  is  much  smaller  than  the  Courant  time  step  at  the 
shock  front.  This  will  occur  when  there  are  regions  of  the  problem  which  require  a 
more  stringent  tim.e  step  than  that  carried  at  the  shock  front.  This  often  happens  in 
a  calculation  which  carries  both  the  fireball  and  the  stress  wave  in  the  ground. 

This  has  led  us  to  consider  the  inclusion  of  additional  explicit  damping.  Some 
care  is  required  because  of  the  potential  for  a  destructive  interaction  with  existing 
damping  mechanisms.  To  prevent  overdamping,  certain  terms  in  the  finite  difference 
equations  were  recentered  to  eliminate  damping  which  was  present  for  time  steps  near 
the  Courant  limit.  There  was  no  modification  of  the  FRAM  scheme  which  provides 
adequate  damping  when  certain  monotonicity  conditions  are  not  satisfied. 

We  then  considered  various  combinations  of  linear  and  quadratic  artificial 
viscosities.  This  exercise  showed  that  linear  artificial  viscosity  was  required  to  control 
the  numerical  pathologies,  but  that  quadratic  viscosity  had  either  small  or  negative 
effects  on  the  solution.  It  appears  that  FRAM  is  adequately  dealing  with  the  problems 
that  would  otherwise  require  quadratic  viscosity,  and  that  additional  quadratic  viscosity 
provides  too  much  damping.  However,  it  appears  that  linear  artificial  viscosity  is 
complementary  to  FRAM,  and  provides  required  damping  where  FRAM  is  Inadequate. 
This  is  believed  to  be  in  the  very  stiff  regime  of  the  equation  state,  where  the 
pressure  is  comparable  or  small  compared  to  the  bulk  modulus. 


1.  Chapman,  M..  ”  FRAM-Nonlinear  Damping  Algorithms  for  the  Continuity 
Equation,"  J.  Comp.  Phys.  44.  84-103,  1981. 
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The  form  of  the  viscosity  used  is  as  follovi/s: 

Q  =  -  p  C,  C  A-v  ,  for  all  A-v 

where 

p  =  material  density 

=  damping  coefficient 

C__  =  local  sound  speed 

A-v  =  velocity  difference 

We  tried  various  values  of  the  damping  coefficient  (C|^),  and  have  found  that  a 
value  of  0.5  is  the  best  compromise  between  controlling  the  damping  and  minimal 
smearing  and  spreading  of  the  wave. 

Appendices  A-D  which  follow  compare  four  calculations  which  clearly 
demonstrate  the  merit  of  explicit  linear  artificial  viscosity  in  addition  to  FRAM.  The 
calculations  model  1-D  spherical  explosions  in  MISTY  ECHO  grout  with  a  starting 
pressure  of  100  Mbar.  The  test  matrix  consists  of  coarse  and  fine  Eulerian  grids,  with 
and  without  the  code  modification.  The  coarse  grid  simulates  the  toning  actually  used 
for  the  Misty  Echo  calculation,  while  the  fine  grid  solutions  are  essentially  fully 
converged. 

Appendix  A  Standard  STREAK  Coarse  zoning 

Appendix  B  Standard  STREAK  Fine  zoning 

Appendix  C  STREAK  with  C  =  0.5  Fine  zoning 

Appendix  D  STREAK  with  C  =  0.5  Coarse  zoning 

The  coarse  solution  without  the  code  modification  produces  waveforms  which 
are  quite  noisy,  and  clearly  questionable.  The  corresponding  fine  solution  is  clearly 
better  behaved,  with  only  a  minor  artifact  at  the  wavefront.  It  is  generally 
unambiguous  how  to  interpret  the  fine  zoned  calculation. 
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With  the  modified  code,  the  fine  zoned  calculation  with  the  modified  code  is 
essentially  identical  to  the  fine  zoned  calculation  with  the  unmodified  code  with  the 
exception  that  overshoots  at  the  wavefront  are  well  controlled.  The  coarse  calculation 
with  the  modified  code  gives  quite  good  agreement  with  the  fine  zoned  calculations, 
with  only  a  minor  overshoot  at  the  wavefront.  The  improvement  over  the  coarse 
calculation  with  the  unmodified  code  is  dramatic. 

In  spite  of  these  results,  linear  artificial  viscosity  must  be  used  with  caution. 
There  are  reports  of  serious  solution  errors  traceable  to  linear  artificial  viscosity.  Two 
such  cases  that  other  at  S-CUBED  have  had  direct  experience  with  are  with  the 
propagation  of  weak  gasdynamic  signals,  and  the  propagation  of  linear  elastic  signals 
in  Lagrangian  codes.  It  is  possible  that  these  difficulties  are  unique  to  Lagrangian 
codes. 


To  test  the  gasdynamic  question,  we  have  performed  a  similar  set  of  test 
calculations  for  a  shock  propagating  in  real  air.  For  the  high  pressure  regime  (shock 
pressures  greater  than  a  few  hundred  bars),  the  linear  artificial  viscosity  is  still  an 
improvement,  albeit  a  small  one.  For  the  weak  shock  case  (shock  pressures  down  to 
1.2  bars),  the  viscosity  is  a  larger  improvement  and  clearly  a  good  feature  to  add  to 
the  code.  The  following  appendices  show  the  results. 


Appendix  E 
Appendix  F 
Appendix  G 
Appendix  H 
Appendix  I 
Appendix  J 


Standard  STREAK  High 
STREAK  with  C  ^  0.5  High 
Standard  STREAK  Low 
Standard  STREAK  Low 
STREAK  with  C  =  0.5  Low 
STREAK  with  C  =  0.5  Low 


pressure 

Fine 

zoning 

pressure 

Fine 

zoning 

pressure 

Coarse  zoning 

pressure 

Fine 

zoning 

pressure 

Fine 

zoning 

pressure 

Coarse  zoning 

In  all  these  cases,  the  version  with  an  artificial  viscosity  clearly  imp'^oves  the 
solution,  and  brings  it  close  to  what  the  standard,  fine  zoned  case  would  suggest 
when  the  obvious  over-ring  is  removed. 
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Even  with  these  .-suits,  however,  we  think  that  the  modified  code  should  be 
used  with  caution  until  further  studies  can  be  completed.  We  need  to  run  a  two- 
dimensior.al  test  and  this  will  be  a  rerun  of  part  of  the  MISTY  ECHO  calculation  to 
see  if  the  problems  clean  up  and  be  sure  there  are  no  unforeseen  problems. 

There  remains  the  problem  of  filtering  numcrica'  oscillations  out  of  existing 
solutions.  The  biggest  difficulty  here  is  estimating  the  true  peak  when  there  has  been 
a  numerical  overshoot.  We  have  taken  two  independent  approaches  to  this  problem, 
and  have  found  the  results  to  be  consistent. 

The  first  approach  is  based  on  the  fact  that  the  observed  peaks  (pressure,  density 
and  velocity)  satisfy  the  equation  of  state  Rankine-Hugoniot  relations,  but  with  a 
shock  speed  which  is  systematically  different  than  the  wave  speed  in  ti.e  calculation. 
If  we  postulate  that  the  calculational  wave  speed  is  correct  and  infer  peak  values  from 
the  equation  of  state  Rankine-Hugoniot  relations,  these  peak  values  a  '  ear  'easonable 
in  every  case  we  examined.  (The  numerical  v^ave  speed  is  inherently  noisy,  and 
required  smoothing.) 

Figures  1  and  2  show  the  correction  factors  as  a  function  of  raw  pressure  and 
velocity  for  the  S-CUBED  MISTY  ECHO  final  prediction. 

As  a  check  on  this  question,  we  ran  a  spherical  ID  calculation  with  similar 
zoning  and  time  steps  as  a  proxy  lor  the  2D  problem.  This  calculation  produ'~ed 
underdamped  waveforms  similar  to  those  in  the  2D  calculations.  We  then  performed 
a  convergence  analysis  by  refining  the  zones  in  the  ID  problem.  The  results  agree 
v^ell  with  the  corrected  data  as  obtained  by  the  procedure  outlined  above. 


5 


i).MtSS\).l(|  0[KJ;)  SA  (d 


) 


d  spoo/d  patJTPOW 


6 


Figure  1  Corrected  pre‘isiire/(  ode  prchsiire  l.xior  versus  rode  irossure  loi  lUJDDY  3 
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Code  Vclocily  \^ca\\ / see) 

Figure  2.  Corrtcled  veloc  ily/(  ode  velocity  (.ictot  versus  code  velocity  lor  lUJODY  J 
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APPENDIX  A 
STANDARD  STREAK 
COARSE  ZONING 
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STREAK  WITH  =  0.5 
LOW  PRESSURE.  COARSE  ZONING 
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